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Using  density  functional  theory  calculations,  many  researchers  have  predicted  that  various  tungsten  nitride 
compounds  Ni_xWx  (x  <  |)  will  be  “ultraincompressible”  or  “superhard”  i.e.,  as  hard  as  or  harder  than 
diamond.  Necessary  conditions  for  such  compounds  are  that  they  have  large  bulk  and  shear  moduli,  greater  than 
approximately  200  GPa,  and  are  elastically  and  vibrationally  stable.  Compounds  with  such  desirable  properties 
also  must  be  energetically  stable  against  decomposition  into  other  compounds.  This  test  for  stability  can  only  be 
found  after  the  determination  of  the  convex  hull  for  Ni_xWx,  which  connects  the  lowest  enthalpy  structures  as 
a  function  of  composition.  Unfortunately,  the  experimental  phase  diagram  of  the  N-W  structure  is  uncertain,  as 
it  is  difficult  to  break  the  N2  bond  to  form  compounds  with  tungsten.  Experiment  also  indicates  that  there  are 
a  large  number  of  partially  filled  sites  in  most  N-W  structures.  This  introduces  computational  difficulties  since 
we  cannot  easily  model  randomly  placed  vacancies.  In  addition,  van  der  Waals  forces  play  a  significant  role  in 
determining  the  structure  of  solid  N2  and  the  nitrogen-rich  compounds.  This  makes  it  difficult  to  determine  the 
relative  energies  of  these  compounds,  as  there  is  no  universally  accepted  density  functional  incorporating  van 
der  Waals  interactions.  The  exact  shape  and  even  composition  of  the  convex  hull  is  dependent  upon  the  choice 
of  density  functional,  even  if  we  only  chose  between  the  local  density  approximation  and  a  generalized  gradient 
functional.  Despite  these  difficulties,  computations  can  determine  much  about  the  ground-state  form  of  the  convex 
hull.  Here,  we  use  high-throughput  calculations  to  map  out  the  hull  and  other  low-energy  structures  for  the  N-W 
system.  The  lowest-energy  structures  all  have  vacancies,  on  the  tungsten  sites  in  hexagonal-based  compounds, 
and  on  both  the  nitrogen  and  tungsten  sites  in  cubic  compounds.  We  find  that  most  of  the  N-W  structures 
proposed  in  the  literature,  both  theoretical  and  experimental,  are  above  the  convex  hull,  in  some  cases  by  over  0.2 
eV/ atom.  One  of  the  ground-state  phases,  N-W  in  the  NbO  structure,  has  relatively  large  bulk  (>300  GPa)  and 
(>200  GPa)  shear  moduli,  and  so  is  a  candidate  superhard  material.  This  will  require  further  investigation. 
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I.  INTRODUCTION 

One  of  the  major  goals  of  electronic- structure  calculations 
is  the  prediction  of  crystal  structures  as  a  function  of 
composition  [1-3].  Determining  the  possible  configurations 
of  a  compound  as  a  function  of  composition  is  the  first  step 
in  determining  its  material  properties  at  equilibrium.  It  is 
also  likely  that  any  pressure-  or  temperature-driven  phase 
transitions  will  be  from  the  equilibrium  structure  to  another 
structure  which  is  close  to  it  in  energy  and  composition.  Such 
calculations  are  of  particular  interest  when  there  is  little  known 
about  the  system  theoretically. 

There  are  a  variety  of  mechanisms  for  this:  searching 
over  a  wide  range  of  known  [4]  and  likely  [5]  structures 
for  the  material  in  question,  searches  starting  from  randomly 
positioned  atoms  [6],  and  even  structures  predicted  from 
apparently  “out  of  the  blue”  [7].  In  the  end  these  techniques 
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produce  a  set  of  metastable  structures,  all  of  which  have  zero 
force  on  the  atoms  in  the  crystal,  zero  stress,  and  no  imaginary 
phonon  modes.  Some  structures  will  be  stable,  that  is,  it  is  not 
energetically  favorable  for  them  to  decompose  into  into  other 
structures. 

When  searching  for  stable  structures,  this  last  point  is  key. 
We  can  find  the  lowest-energy  structure  of,  e.g.,  composition 
A2B5,  but  we  cannot  tell  if  it  is  truly  stable  until  we  show  that 
it  cannot  decompose  into  2A  +  5B,  or  A2B2  +  3B,  or  some 
other  combination  of  structures.  In  other  words,  to  completely 
determine  structural  stability  of  a  compound,  we  must  search 
over  all  possible  compositions. 

Once  this  is  done  there  is  still  the  question  of  whether 
or  not  we  have  done  the  appropriate  calculation.  Most 
electronic- structure  calculations  are  done  using  the  Kohn- 
Sham  ansatz  [8]  to  density  functional  theory  (DFT)  [9].  Since 
the  exact  density  functional  is  not  known,  we  necessarily 
use  approximate  forms.  It  is  not  a  priori  certain  that  the 
local  density  approximation  (LDA)  [10,11]  will  give  the 
same  results  as  a  generalized-gradient  functional  such  as 
Perdew-Burke-Ernzerhof  (PBE)  [12].  If  van  der  Waals  forces 
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are  important,  we  should  use  some  variety  of  van  der  Waals- 
enabled  functional  [13].  If  multiple  functionals  for  a  given 
system  predict  the  same  stable  structures  and  compositions, 
then  it  is  likely  that  these  compounds  will  form.  If  different 
functionals  produce  different  results,  then  we  cannot  decide 
the  issue  on  the  basis  of  current  density  functional  theory  and 
must  appeal  to  more  accurate  techniques  or  to  experiment. 

As  an  example,  we  look  at  tungsten  nitrogen  system.  There 
are  many  experimental  studies  showing  evidence  of  stable 
Ni_xWx  compounds  [14-21],  and  there  is  no  assessed  phase 
diagram  diagram  available  [22].  Indeed,  until  recently  [21]  it 
was  very  difficult  to  make  nitrogen-rich  Ni_xWx  compounds 
since  tungsten  does  not  dissolve  appreciable  amounts  of  N2 
[23].  In  addition,  above  600  °C  the  compound  usually  referred 
to  as  tungsten  nitride  decomposes  to  NW2,  which  itself  loses 
nitrogen  above  1000  °C  [23]. 

Tungsten  carbide  is  of  interest  as  a  structural  material,  or  at 
least  a  coating,  because  of  computational  studies  which  claim 
that  some  compounds  Ni_xWx  with  v  <  ^  are  “superhard”, 
or  at  least  have  very  stiff  elastic  constants  and  so  are  superhard 
candidates  [24-30]. 

While  all  these  calculations  show  that  the  favored  phase 
of  Ni_xWx  is  stable  against  decomposition  into  N2  and 
tungsten,  none  provide  a  detailed  calculation  of  the  possible 
stable  and  metastable  structures  of  N\-XWX  as  a  function  of 
tungsten  concentration  v.  It  is  thus  difficult  to  determine  if 
these  structures  will  form  experimentally,  or  to  determine  if 
structures  which  do  form  are  stable  or  only  metastable. 

Our  original  motivation  for  this  study  was  the  possibility 
that  N4W  (or  N0.8W0.2)  might  be  a  superhard  metal  [30] 
and  can  be  in  fact  easily  fabricated.  As  noted  above,  this 
requires  knowledge  of  the  entire  equilibrium  phase  diagram 
for  N 1  _x  Wx .  We  did  this  using  AFLOW  [3 1  ] ,  a  high-throughput 
front  end  for  electronic- structure  calculations  [1].  AFLOW 
allows  us  to  quickly  examine  a  selected  range  of  structures 
using  high-performance  supercomputers  and  modem  den¬ 
sity  functional  electronic- structure  techniques.  The  AFLOW 
prototype  database,  which  was  originally  used  to  describe 
intermetallic  alloys  [32-36],  can  easily  be  enlarged  [37].  We 
were  therefore  able  to  include  ionic  and  covalent  structures 
for  systems  similar  to  N-W.  These  include  borides,  carbides, 
oxides,  and  other  nitrides.  We  also  invented  many  structures 
to  mimic  the  random  pattern  of  vacancies  on  both  the  tungsten 
and  nitrogen  sites  that  have  been  observed  under  experimental 
conditions. 

AFLOW  allows  the  use  of  a  variety  of  different  density 
functionals.  We  did  our  primary  calculations  using  the  Perdew- 
Burke-Ernzerhof  (PBE)  generalized-gradient  functional  [12]. 
For  low-energy  systems  we  also  used  the  local  density  approx¬ 
imation  [10,11]  (LDA).  We  also  approximated  van  der  Waals 
forces,  which  are  important  in  nitrogen-rich  structures,  using 
the  vdW-DF2  functional  [38].  As  we  shall  see,  predictions  for 
ground-state  structures  are  strongly  dependent  on  the  choice 
of  functional. 

Here,  we  are  searching  for  structures  of  composition 
Ni_xWx  which  are  stable,  or  are  metastable  with  a  reasonable 
probability  of  being  found  experimentally.  A  metastable 
compound  must  be  vibrationally  stable,  i.e.,  have  no  imaginary 
phonon  frequencies  as  well  as  elastic  constants  which  satisfy 
the  Born  criteria  [39].  Tme  stability  requires  an  additional 
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step:  consider  the  decomposition  of  a  compound  N^Wg  into 
the  ground-state  forms  of  its  end-point  components,  aN2  [40], 
and  body-centered-cubic  tungsten  [41].  This  decomposition 
requires  a  change  in  enthalpy 

AH(NaWb)  =  [E( NaW5)  -  B  E( W) 

— (1/2)A  E(aN2)\/(A  +  B)  <  0,  (1) 

where  A//(NaW5)  is  the  formation  energy  per  atom 
of  a  compound  with  this  stoichiometry,  E(NaWb)  is  the 
energy/formula  unit  of  the  compound,  E( W)  is  the  equilibrium 
energy  per  atom  of  body-centered-cubic  tungsten,  and  E(otN2) 
is  the  equilibrium  energy  per  molecule  of  solid  aN2.  A  H  will 
be  negative  if  it  is  energetically  possible  for  the  compound  to 
be  formed  spontaneously  from  solid  N2  and  W.  The  compound 
will  be  completely  stable  if  it  is  not  energetically  favorable  for 
it  to  decompose  into  any  other  N-W  compound.  This  requires 
that  the  compound  be  on  the  convex  hull  determined  from  the 
plot  of  AH  versus  tungsten  concentration  v  =  N/(M  +  N ) 
over  all  possible  compounds  in  the  system  [42]. 

The  purpose  of  the  paper  is  to  determine  AH  for  a 
sufficient  number  of  N-W  compounds  so  that  we  have  a 
good  approximation  to  the  shape  of  the  convex  hull  and  the 
structures  forming  it,  and  to  determine  the  electronic,  elastic, 
and  vibrational  properties  of  compounds  on  or  near  the  convex 
hull,  looking  for  superhard  candidates.  An  estimate  of  the 
actual  hardness  of  these  materials  is  beyond  the  scope  of  this 
paper,  but  in  general  will  follow  the  procedures  in  Ref.  [30]. 

The  paper  is  organized  as  follows.  Experimental  and  previ¬ 
ous  theoretical  work  is  shown  in  Sec.  II.  Our  computations  are 
described  in  Sec.  III.  Our  results  are  shown  in  Sec.  IV,  which 
is  broken  up  into  several  subsections,  including  /3  (cubic) 
and  8  (layered  hexagonal)  compounds,  compounds  containing 
N2  dimers,  compounds  resembling  Si02,  and  miscellaneous 
structures.  Finally,  we  discuss  our  results  in  Sec.  V. 

Our  study  of  tungsten  nitride  required  searching  over  a 
large  number  of  known  experimental  structures  [43-47],  and 
even  inventing  structures  which  had  vacancy  patterns  related 
to  those  seen  experimentally.  The  Supplemental  Material  [48] 
contains  crystallographic  information  for  108  of  the  lowest- 
energy  structures,  along  with  additional  computational  details. 
In  all,  we  studied  over  500  structures  using  the  high-throughput 
AFLOW  platform.  Structural  parameters  for  these  additional 
structures  are  available  on  request. 

We  always  refer  to  the  composition  of  tungsten  nitride  in 
the  form  Ni_xWx,  where  v  is  the  atomic  fraction  of  tungsten, 
or,  for  stoichiometric  compounds,  N^W#,  where  A  and  B 
are  integers.  Various  other  authors  reverse  the  order  of  these 
compounds,  e.g.,  “r-W2N3”  [21].  To  avoid  confusion,  we  place 
these  designations  in  quotes.  In  the  Supplemental  Material 
[48],  all  structures  are  ordered  by  tungsten  concentration  v. 

II.  EXPERIMENTAL  AND  THEORETICAL  BACKGROUND 

Wriedt  [22]  summarizes  the  experimental  N-W  data  known 
before  1990.  The  cubic  phase  is  essentially  the  face- 
centered-cubic  NaCl  structure  (structure  #67  in  the  Supple¬ 
mental  Material  [48]),  with  vacancies  on  the  nitrogen  site 
leading  to  compounds  of  the  form  Ni_xWx  with  x  >  ^.  The 
cubic  lattice  constant  for  all  of  these  phases  is  approximately 
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4.13  A.  These  structures  were  originally  studied  by  Hagg 
[14],  and  later  by  Khitrova  and  Pinsker  [17-19].  Kiessling 
and  co-workers  [15,16]  found  a  similar  phase  which  they 
designated  y.  This  structure  also  seems  to  be  related  to  the 
NaCl  structure,  but  with  x  ~  1. 

The  hexagonal  8  phase  structures  were  extensively  studied 
by  Khitrova  and  Pinsker  [18,19].  They  defined  six  distinct 
crystallographic  structures  8^IV  and  8\,IV ,  which  we  will 
discuss  in  more  detail  in  the  following.  Compositions  Ni_xWx 
are  said  to  range  from  x  =  0.33  to  0.67.  As  NaCl  is  the  base 
structure  for  the  ft  phase,  we  can  consider  tungsten  carbide 
[49]  (structure  #61  in  the  Supplemental  Material  [48])  as 
the  basis  for  most  of  the  8  phase  structures.  These  structures  are 
composed  of  stacked  triangular  planes  of  tungsten  or  nitrogen 
atoms.  Unlike  the  ft  phases,  vacancies  appear  only  in  the 
tungsten  planes. 

In  addition  to  bulk  structures,  Shen  and  Mai  [50]  have 
grown  thin-film  N-W  films.  They  found  some  evidence  of  the 
ft  phase  at  NW2,  but  most  of  the  films  were  amorphous. 

These  compounds  are  difficult  to  make  and  to  stabilize.  As 
noted  by  Toth  [23],  “tungsten  does  not  dissolve  appreciable 
amounts  of  nitrogen  as  a  terminal  solid  solution,”  and  must 
be  formed  at  temperatures  under  800  °C.  Recently,  however, 
Wang,  Yu,  Lin  et  al.  [21]  (hereafter  referred  to  as  WYL+)  were 
able  to  use  solid-state  ion  exchange  and  nitrogen  degassing 
under  pressure  to  produce  several  tungsten  nitride  phases  with 
compositions  ranging  from  N3W2  to  N2W3  (0.4  <  x  ^  0.6). 
Although  they  do  not  state  this,  all  of  these  structures  can 
be  classified  as  either  fi-NW  or  <5-NW,  as  defined  above.  We 
will  have  more  to  say  about  these  structures  as  we  discuss  our 
computational  results. 

Computational  studies  of  the  tungsten-nitride  system  fol¬ 
lowed  in  the  wake  of  the  development  of  high-speed  computers 
and  computationally  efficient  density  functional  codes.  We  will 
discuss  some  of  these  papers  in  more  detail  in  the  next  section. 
For  now,  we  briefly  describe  their  results.  These  researchers 
studied  specific  crystal  structures  of  tungsten  nitride,  but  none 
have  examined  structures  over  a  large  range  of  compositions. 

Kroll,  Schroter,  and  Peters  [24]  considered  the  nickel 
arsenide  structure  [51]  (#59  in  the  Supplemental  Material  [48]) 
to  be  the  ground  state  of  NW.  They  also  looked  at  several 
possible  structures  for  N2W,  including  baddeleyite  [52]  (#20), 
which  they  found  to  have  the  lowest  energy,  brookite  [53] 
(#18),  and  cotunnite  [54]  (#23),  which  they  found  to  be  a 
high-pressure  (>35  GPa)  phase  of  N2W. 

Suetin,  Shein,  and  Ivanovski  [25]  modeled  tungsten  nitride 
by  assuming  it  to  be  in  either  the  tungsten  carbide  (WC,  #61  in 
the  Supplemental  Material  [48])  or  the  sodium  chloride  (NaCl, 
#67)  structure.  They  computed  elastic  constants,  isotropic  bulk 
moduli  in  the  Reuss  and  Voigt  approximations  [55],  and  the 
corresponding  isotropic  Young’s  modulus  and  Poisson’s  ratio. 
The  NaCl  structure  is  elastically  unstable,  with  C44  <  0.  The 
WC  compound,  on  the  other  hand,  was  found  to  be  rather  stiff, 
with  a  shear  modulus  of  about  150  GPa.  The  enthalpy  AH 
was  not  calculated,  so  there  could  be  no  estimate  of  stability 
against  dissociation  into  N2  and  bcc  tungsten. 

Wang,  Li,  Li,  Xu  et  al.  [26]  (hereafter  WLLX+)  used  an 
evolutionary  method  to  find  two  related  N2W  structures  (#11 
and  #12)  which  had  shear  moduli  >200  GPa.  These  structures, 
which  consisted  of  vertically  aligned  N2  dimers  alternating 
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with  hexagonal  tungsten  planes,  have  AH  <  0,  as  well  as 
elastic  constants  satisfying  the  Bom  criteria  [39].  In  addition, 
both  structures  are  insulators,  with  a  band  gap  on  the  order  of 
1  eV. 

Benhai,  Chunlei,  Xuanyu,  Qiuju,  and  Dong  [27]  looked  at 
the  elastic  behavior  of  a  variety  of  NW  (x  =  \)  structures,  at 
pressures  up  to  100  GPa.  They  found  that  the  NiAs  structure 
was  preferred  over  the  WC  at  all  pressures.  They  did  not 
calculate  the  change  in  enthalpy  relative  to  the  N2  and  bcc 
W  end  points,  nor  did  they  consider  some  of  the  lower-energy 
structures  we  report  in  the  following. 

Song  and  Wang  [29]  studied  NW  in  the  WC  structure,  N2W 
in  the  CoSb2  structure  (#20,  although  they  refer  to  it  as  IrP2), 
and  N3W  in  the  P3TC  structure  [56]  (#10),  using  the  LDA 
functional.  The  found  a  negative  value  of  AH  for  all  three 
compounds,  although  it  is  not  clear  what  reference  point  they 
used  for  pure  nitrogen.  In  all  cases,  they  found  shear  moduli 
below  200  GPa  but  bulk  moduli  greater  than  300  GPa. 

Suetin,  Shein,  and  Ivanovski  [57]  modeled  N2W  in  what  we 
believe  to  be  the  CTi2  stmcture  [58],  a  cubic  supercell  of  the 
NaCl  structure  with  ordered  vacancies  on  one  of  the  sublattices 
(see  structure  #95  in  the  Supplemental  Material  [48]).  This  is 
an  approximation  to  the  cubic  NW2  phase  found  by  Hagg  [14] . 
Although  they  compared  the  electronic  structure  to  that  found 
for  the  WC  and  NaCl  structures  in  their  earlier  paper  [25],  they 
give  no  information  about  the  stability  of  this  structure. 

Du,  Wang,  and  Lo  [28]  looked  at  a  tetragonal  analog  (#26) 
to  the  hexagonal  structures  found  by  WLLX+.  They  found 
it  to  be  stable  compared  to  these  structures  above  150  GPa, 
but  never  stable  compared  to  the  cotunnite  structure  (#23).  Li, 
Zhai,  Fu  et  al  [59]  computed  the  pressure  dependence  of  the 
elastic  constants  of  these  structures,  and  found  that  both  are 
elastically  stable  up  to  pressures  of  200  GPa. 

In  addition  to  their  experimental  work,  WYL-f-  [21]  also 
used  first-principles  calculations  to  determine  the  elastic 
constants  of  the  three  stoichiometric  structures  they  studied,  as 
well  as  what  they  call  the  c-BN  structure  (actually  zinc-blende 
[60-62]  #65).  Bulk  moduli  all  exceed  350  GPa,  and  shear 
moduli  are  estimated  at  around  180  GPa,  but  reach  390  GPa 
for  the  c-BN  structure.  They  did  not  otherwise  discuss  the 
possible  stability  of  these  phases. 

Ay  din,  Ciftci,  and  Tatar  [30]  (hereafter  ACT)  modeled  N4W 
using  the  starting  point  of  the  ReP4  structure  [63]  (#7).  They 
found  this  structure  to  be  mechanically  stable,  and  AH  < 
0.  The  shear  modulus  of  this  structure  was  approximately 
200  GPa,  and  they  estimated  the  hardness  of  the  material  to  be 
on  the  order  of  c-BN.  If  it  could  be  made,  then,  it  would  most 
likely  be  a  superhard  material.  We  will  discuss  our  findings  for 
this  structure  in  detail  in  the  following. 

Zhang,  Yan,  Wei,  and  Wang  [64]  calculated  the  equilibrium 
lattice  parameter  and  elastic  constants  of  WYL+’s  “c- W3N4” 
structure  [21],  which  has  the  prototype  [65]  S3U4  (#47).  Later 
Liu,  Wang,  Zhou,  and  Chang  [66]  looked  at  the  pressure 
dependence  of  the  elastic  constants  and  phonon  frequencies 
of  the  same  structure.  This  structure  has  very  large  bulk  and 
shear  moduli,  but  we  shall  see  that  it  is  far  above  the  convex 
hull  of  the  tungsten  nitride  phase  diagram. 

Liu,  Zhou,  Gall,  and  Khare  [67]  recently  computed  the 
elastic  constants  of  transition-metal  nitrides  in  the  NbO 
structure  [68].  Their  calculation  for  NW  showed  it  to  be  very 
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stiff,  with  its  bulk  modulus  over  300  GPa  and  the  shear  modulus 
above  200  GPa.  We  will  discuss  this  structure  further  below. 

Summarizing  the  experimental  and  theoretical  results,  we 
find  that  the  low-lying  structures  of  tungsten  nitride  fall  into 
one  of  five  classes,  the  first  two  roughly  corresponding  to 
structures  discussed  by  Schonberg  [49]: 

(1)  phases.  These  structures  are  supercells  of  the  cubic 
NaCl  structure  (#67  in  the  Supplementary  Material  [48]). 
Removing  atoms  from  the  supercells  generally  leads  to  a 
lower-energy  cell.  Examples  are  the  S3U4  structures  (#47 
and  #78),  which  have  either  one  N  or  one  W  vacancy  in  an 
eight-atom  supercell,  the  NbO  [68]  (#53)  structure,  which  has 
second-neighbor  N  and  W  vacancies,  and  the  CTi2  structure 
[58]  (#95),  which  approximates  the  NW2  ft  phase  found  by 
Hagg  [14].  For  our  purposes  we  will  treat  the  similar  y-NW 
structure  [15,16]  as  part  of  the  ft  phase. 

(2)  8  phases.  In  general  these  are  hexagonal  or  trigonal 
unit  cells  containing  alternating  triangular  planes  of  nitrogen 
and  tungsten.  The  tungsten  planes  generally  contain  vacancies 
[18,21].  The  experimental  composition  of  these  phases  ranges 
from  about  [21]  N4W3  to  NW2  [18]. 

(3)  N2  phases.  In  these  structures  the  N2  molecules  stay 
intact,  possibly  within  a  tungsten  matrix.  The  only  exper¬ 
imental  examples  we  have  of  these  phases  are  the  several 
varieties  of  solid  N2  [40] .  Computations  have  found  structures 
with  composition  N2W  (Ref.  [26],  #11  and  #12)  and  N4W 
(Ref.  [30],  #7)  that  appear  to  be  at  least  metastable. 

(4)  SiC>2  phases.  These  have  not  been  seen  experimentally, 
and  are  only  exceptionally  low  in  energy  when  using  the  vdW- 
DF2  van  der  Waals  functional.  In  these  structures,  the  nitrogen 
atoms  form  a  tetrahedron  around  each  tungsten  atom,  and  each 
nitrogen  bonds  with  two  tungsten  atoms.  A  number  of  such 
structures  appear  in  the  literature  [69].  Here,  we  examine  a 
few  of  the  more  common  ones. 

(5)  Miscellaneous  structures.  This  category  includes  those 
structures  which  do  not  fit  into  one  of  the  categories  defined 
above,  including  a  derivative  of  the  FeB4  compound  proposed 
by  Van  der  Geest  and  Kolmogorov  [70]  (#8),  and  the 
surprisingly  low-energy  MoS2  structure  (#89).  None  of  these 
have  been  seen  experimentally. 

III.  METHODS 

We  began  our  search  for  the  convex  hull  of  the  N-W  system 
by  using  AFLOW  [1,31,71]  to  quickly  and  efficiently  search 
through  a  large  database  of  structures.  As  the  original  AFLOW 
prototype  database  was  for  binary  metallic  alloys  [32-36], 
we  have  extended  it  to  include  over  100  new  structures. 
These  include  nitrides,  oxides,  borides,  and  carbides,  as  well 
as  supercells  of  standard  structures  with  atoms  removed  to 
mimic  the  random  patter  of  vacancies  seen  experimentally. 
The  crystallographic  information  for  the  important  structures, 
including  all  of  those  discussed  in  Secs.  II  and  IV  are  described 
in  the  Supplemental  Material  [48]. 

Electronic- structure  calculations  were  done  using  the  Vi¬ 
enna  ab  initio  simulation  package  (VASP)  [72,73],  including 
core- state  effects  via  the  VASP  implementation  [74]  of  the 
projector  augmented- wave  (PAW)  method  [75]. 

aflow’s  default  is  to  use  the  Perdew-Burke-Ernzerhof 
(PBE)  implementation  [12]  of  the  generalized-gradient 
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approximation  (GGA)  to  density  functional  theory  (DFT) 
[8,9].  There  are,  however,  substantial  differences,  especially 
in  the  estimates  of  enthalpy  changes,  between  PBE  results  and 
those  from  the  LDA  [10,11].  This  was  shown  by  WLLX+ 
[26],  and  we  shall  see  further  examples  in  the  following. 
Accordingly,  we  calculated  the  convex  hull  using  both  LDA 
and  PBE  functionals  within  VASP. 

Molecular  nitrogen,  ground  state  [40]  aN2  (#1  or  #2  in  the 
Supplemental  Material  [48]),  is  a  van  der  Waals  solid  [76], 
and  neither  the  LDA  nor  the  PBE-GGA  correctly  predict  its 
lattice  constant.  We  studied  the  effect  of  van  der  Waals  forces 
by  utilizing  the  VASP  implementation  [13,77]  of  the  vdW-DF2 
functional  [38]  This  functional  was  developed  to  handle  the 
van  der  Waals  interaction  in  general  geometries,  so  we  apply 
it  over  the  entire  range  of  compositions,  calculating  a  third 
convex  hull  and  phase  diagram. 

We  used  the  VASP-supplied  LDA  and  PBE  PAW  potentials 
for  nitrogen  and  tungsten  (specifically,  May  2000  LDA  and 
April  2002  PBE  PAW  potentials  for  N,  and  the  July  1998  LDA 
and  September  2000  tungsten  W_pv  PAWs).  All  calculations 
use  a  kinetic  energy  cutoff  of  560  eV,  which  is  40%  larger  than 
the  suggested  cutoff  for  nitrogen  (400  eV).  Unless  otherwise 
stated,  we  let  the  AFLOW  package  determine  the  k -point  mesh 
for  each  structure.  This  is  usually  set  to  use  approximately  the 
same  density  of  k  points  in  reciprocal  space  for  all  structures. 

Since  there  are  a  wide  variety  of  structures  in  this  system, 
both  metallic  and  insulating,  we  made  sure  that  the  k-point 
mesh  was  dense  enough  to  determine  the  electronic-structure 
energy  to  better  than  1  meV / atom.  For  the  NaCl  structure,  e.g., 
a  15  x  15  x  15  Monkhorst-Pack  grid  [78]  was  originally  used. 
For  Wang  et  al. ’ s  P  6m  2  N2W  structure,  we  used  a  15  x  15  x  9 
mesh. 

In  the  initial  VASP  calculations  performed  via  AFLOW,  all 
structures  are  considered  to  be  spin  polarized,  but  the  starting 
magnetization  vanished  as  self-consistency  was  reached.  The 
final  calculations  were  all  done  assuming  no  moment. 

The  MedeA®  software  system  [79]  was  used  to  drive  VASP 
in  calculations  of  the  phonon  spectra  and  for  some  other 
calculations.  Elastic  constants  were  computed  by  taking  finite 
strains  [80,81]  and  determining  the  slope  of  the  corresponding 
stress-strain  curve  (as  implemented  in  the  MedeA®  package, 
or  using  the  native  VASP  option).  Phonon  frequencies  were 
determined  via  the  frozen-phonon  approximation  using  the 
MedeA®  package. 

IV.  COMPUTATIONAL  RESULTS 

The  surprising  thing  about  our  computational  work  is  that 
so  little  of  it  agrees  with  either  the  experimental  or  theoretical 
work  discussed  above.  This  section  details  our  results,  using 
three  different  density  functionals.  In  the  next  section,  we  will 
discuss  the  disagreement  between  these  results  and  previous 
work. 

We  computed  the  enthalpy  of  formation  AH  [Eq.  (1)]  for 
over  500  Ni_xWx  structures  using  the  AFLOW  mechanism  for 
high-throughput  calculations,  with  VASP  as  our  computational 
engine  [72,73,82].  We  used  the  PBE  density  functional  [12] 
as  our  screening  tool,  but  for  “interesting”  structures,  i.e., 
structures  within  approximately  0.2  eV  of  the  tie  line,  we 
also  looked  for  the  minimum  energy  structure  using  the  LDA 
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FIG.  1.  (Color  online)  The  relative  enthalpy  (1)  of  Ni_xWx  for  the  108  structures  listed  in  the  Supplemental  Material  plotted  versus 
increasing  tungsten  concentration.  These  points  were  calculated  using  the  local  density  approximation  (LDA)  with  VASP.  Red  squares  mark  the 
P -phase  structures,  dark  blue  diamonds  the  8  phases,  purple  circles  the  N2  (nitrogen  dimer)  phases,  gray  triangles  the  Si02  phases,  and  green 
pentagons  the  remaining  miscellaneous  phases.  The  labels  correspond  to  the  numbers  of  the  structures  in  the  Supplemental  Material  [48]. 


[10,11]  and  vdW-DF2  [38]  functionals.  We  also  included 
all  of  the  stoichiometric  structures,  both  experimental  and 
theoretical,  from  Sec.  II.  In  the  case  of  experimental  structures 
with  randomly  ordered  vacancies,  we  constructed  supercells  of 
the  base  structure,  removing  nitrogen  and/or  tungsten  atoms  as 
needed  to  reach  the  desired  composition.  This  resulted  in  many 
high-energy  structures  which  we  will  not  describe  further,  but 
the  procedure  also  found  structures  which  shape  the  convex 
hull  of  the  energy  diagrams. 

The  resulting  phase  diagrams  for  these  108  structures  are 
shown  in  Figs.  1  (LDA  functional),  2  (PBE  functional),  and  3 
(vdW-DF2  functional).  The  diagrams  are  decidedly  different, 
especially  for  tungsten  concentrations  x  <  0.4.  The  shape, 
depth,  and  even  composition  of  the  convex  hull  changes  as  we 
move  from  LDA  to  PBE  to  vdW-DF2  functionals.  Differences 
between  functionals  will  be  discussed  in  the  comments  on  the 
individual  structures,  below. 

In  the  following,  we  describe  in  detail  the  most  interesting 
of  the  108  structures  listed  in  the  Supplemental  Material  [48], 
paying  particular  attention  to  how  our  predictions  change  with 
change  in  the  choice  of  density  functional. 


A.  End  points  of  the  phase  diagram 

First  consider  the  end-point  structures:  The  lowest-energy 
N2  structures  are  described  by  Donohue  [40]  and  Mills, 
Oligner,  and  Cromer  [83],  comprising  structures  #l-#5  in 
the  Supplemental  Material  [48]  tables.  We  ignore  the  higher- 
energy  nonmolecular  structures  of  nitrogen  such  as  cg-N  [84], 
considering  only  solid  N2  crystals,  where  the  molecules  are 


bound  by  their  mutual  van  der  Waals  attraction  [76].  Since 
LDA  and  PBE  do  not  explicitly  include  long-range  van  der 
Waals  forces,  we  might  expect  some  difficulty  in  describing 
the  ground  state. 

Table  I  shows  the  calculation  of  the  equilibrium  structure 
of  aN2  in  the  Pa 3  structure,  which  is  one  of  two  proposed 
structures.  The  equilibrium  lattice  constant  and  bulk  modulus 
were  determined  by  fitting  energy  versus  volume  data  to  a 
third-order  Birch  fit  [86]  and  the  internal  parameter  x  was 
determined  by  minimizing  the  forces  on  the  atoms  at  fixed 
lattice  constant  [87].  While  all  three  functionals  correctly 
determine  the  N2  bond  length  at  1 . 10  A,  none  does  particularly 
well  in  determining  the  equilibrium  lattice  constant  and  bulk 
modulus. 

All  N2  structures  are  essentially  degenerate,  differing  by  no 
more  than  0.03  eV  per  atom  as  we  go  from  a^2  to  6N2.  Thus 
current  versions  of  DFT  and  electronic- structure  algorithms, 
even  including  van  der  Waals  interactions,  cannot  accurately 
determine  the  ground  state  of  N2.  We  will  study  this  in  more 
depth  in  a  future  paper,  but  for  now  we  will  simply  note 
that  DFT  calculations  for  pure  nitrogen,  and  presumably  for 
nitrogen-rich  compounds,  differ  considerably  depending  on 
the  choice  of  functional.  One  consequence  of  this  is  that  the 
minimum  entropy  for  this  system  ranges  from  —0.35  to 
—0.75  ev/atom,  depending  on  the  choice  of 
functional. 

Tungsten  is  naturally  found  in  two  forms  [41],  a— W 
(#106),  the  ground-state  body-centered-cubic  structure,  and 
/3-W  (#107),  commonly  referred  to  by  its  Strukturbericht 
[43]  designation  A 15.  For  this  calculation  we  also  considered 
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FIG.  2.  (Color  online)  The  relative  enthalpy  (1)  of  Ni_xWx  for  the  108  structures  listed  in  the  Supplemental  Material  [48]  plotted  versus 
increasing  tungsten  concentration.  These  points  were  calculated  using  the  Perdew-Burke-Ernzerhof  (PBE)  functional  with  VASP.  Note  that  the 
enthalpy  of  structure  #7  is  above  the  upper  limit  of  the  graph.  The  meaning  of  the  symbols  is  described  in  the  caption  to  Fig.  1. 

the  face-centered-cubic  structure  (#108).  Table  II  shows  difference  predicted  between  the  ground- state  bcc  structure 
the  equilibrium  lattice  constants  for  tungsten  as  a  function  and  the  other  two  structures.  As  usual,  LDA  underestimates  the 
of  structure  and  density  functional,  as  well  as  the  energy  equilibrium  lattice  constant,  while  PBE  overestimates  it,  and 
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FIG.  3.  (Color  online)  The  relative  enthalpy  (1)  of  Ni_xWx  for  the  108  structures  listed  in  the  Supplemental  Material  [48]  plotted  versus 
increasing  tungsten  concentration.  These  points  were  calculated  using  the  van  der  Waals  vdW-DF2  functional  with  VASP.  The  meaning  of  the 
symbols  is  described  in  the  caption  to  Fig.  1. 
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TABLE  I.  Equilibrium  lattice  constant,  internal  parameter,  and 
equilibrium  bulk  modulus  for  the  Pa3  structure  (#1)  of  ofN2, 
as  determined  using  various  density  functionals  and  compared  to 
experiment  [45,85].  The  lattice  is  simple  cubic,  and  the  nitrogen  atoms 
sit  on  the  (8c)  Wyckoff  position,  which  has  one  internal  parameter 
x.  The  equilibrium  lattice  constant  and  bulk  modulus  are  determined 
from  a  fourth-order  Birch  fit.  The  quantity  d( N-N)  is  the  length  of 
the  nitrogen-nitrogen  bond  in  N2  molecules. 


Functional 

LDA 

PBE 

vdW-DF2 

Expt.  [45] 

a  (A) 

5.223 

6.187 

5.511 

5.659 

0.061 

0.052 

0.058 

0.056 

d( N-N)  (A) 

1.10 

1.11 

1.11 

1.10 

K0  (GPa) 

5.70 

0.788 

4.69 

1.2  [85] 

the  vdW-DF2  lattice  constant  is  even  larger.  All  functionals 
give  approximately  the  same  energy  difference  between  the 
three  structures,  indicating  that  the  van  der  Waals  contribution 
to  the  energy  is,  as  one  would  expect,  negligible  for  pure 
tungsten. 

B.  p -phase  structures 

The  /3 -phase  structures  [14]  (including  the  y-NW  phase  de¬ 
scribed  by  Kiessling  and  co-workers  [17-19])  can  be  described 
as  tungsten  nitride  in  the  NaCl  structure  with  vacancies  on 
selected  nitrogen  and/or  tungsten  sites.  In  the  extreme  limits, 
this  leads  to  the  NaCl  structure  (#67  in  the  Supplemental 
Material  [48]),  with  no  vacancies,  and  a  face-centered-cubic 
structure  of  tungsten  (#108),  with  all  nitrogen  sites  empty. 
Experiments  [14-19,21,23]  generally  assume  that  the  tungsten 
sites  in  ft  (and  y)  N-W  are  fully  occupied,  and  that  changes 
in  stoichiometry  are  controlled  by  vacancies  on  the  nitrogen 
sites,  with  possible  tetrahedral  nitrogen  interstitials  for  the 
nitrogen-rich  phases.  As  we  will  see  in  the  following,  density 
functional  calculations  do  not  support  this  point  of  view. 

The  lowest-energy  ft -phase  structure,  independent  of  the 
choice  of  functional,  is  NbO  (#53)  [68],  a  simple  cubic 
supercell  of  NaCl  with  vacancies  at  the  corners  and  in  the 
center  of  the  cube,  as  shown  in  Fig.  4.  In  itself  this  should  not  be 
particularly  surprising,  as  NW  and  NbO  have  the  same  number 
of  electrons  in  the  valence  band.  However,  no  experimental 
papers  claim  discovery  of  the  NbO  structure.  The  closest 
structure  to  NbO  was  described  by  Kiessling  and  Peterson 
[16]  as  a  structure  with  space  group  Pm3m  with  metal  atoms 

TABLE  II.  Equilibrium  lattice  constants  in  A  and  relative  energy 
differences  in  eV/ atom  for  the  body-centered-cubic  (#106),  A15 
(#107),  and  face-centered-cubic  structures  (#108),  compared  to 
experimental  information  where  available  [41]. 
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0.489 

PBE 

3.1893 
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FIG.  4.  (Color  online)  Tungsten  nitride  in  the  NbO  structure  [68] 
(#53  in  the  Supplemental  Material  [48])  constructed  by  taking  equally 
spaced  vacancies  from  both  sites  of  the  NaCl  structure.  The  PBE 
functional  predicts  the  equilibrium  lattice  constant  to  be  4.131  A. 


on  the  (3c)  Wyckoff  sites.  Instead  of  nitrogen  occupying  the 
(3d)  Wyckoff  sites,  as  would  be  found  in  NbO,  they  found  both 
oxygen  and  nitrogen  on  the  (lb)[l/2  1/2  1/2]  and  (3d)  sites, 
equivalent  to  the  S3U4  structure  [65].  WYL+  [21]  described 
a  structure  which  they  call  “C-W3N4”  (x  =  0.429)  containing 
no  oxygen,  which  also  takes  on  the  S3U4  structure.  While  we 
did  examine  this  structure  (#47),  we  did  not  find  it  particularly 
close  to  the  ground-state  hull  using  any  functional,  although  it 
does  have  a  lower  value  of  AH  than  the  NaCl  structure. 

The  NbO  structure  for  tungsten  nitride  was  recently 
examined  by  Liu,  Zhou,  Gall,  and  Khare  [67],  who  computed 
its  equilibrium  lattice  constants  and  elastic  constants.  They  did 
not,  however,  compare  its  stability  to  the  NaCl  structure,  nor 
find  its  place  in  the  overall  tungsten  nitride  phase  diagram. 

Our  results  show  that  many  low-energy  /3-type  structures 
in  the  N-W  system  can  be  constructed  starting  with  supercells 
of  the  NaCl  structure  and  removing  selected  N  and  W  atoms. 
The  general  rule  we  found  is  that  the  resulting  structure  will 
have  lower  energy  than  NaCl  if  the  resulting  vacancies  are 
not  nearest  neighbors  (separation  a/\fl  for  like  atoms,  all 
for  unlike  atoms).  This  includes  structures  with  large  atomic 
relaxations,  e.g.,  the  M02N  structures  [90]  (#19  and  #92),  and 
structures  which  only  allow  relaxation  of  the  lattice  constant, 
e.g.,  S3U4  (#47  or  #78)  and  the  ground-state  NbO  structure. 

The  elastic  constants  of  the  NbO  phase  are  given  in 
Table  III,  along  with  the  results  of  Liu  et  al.  [67].  Each 
functional  predicts  this  structure  to  be  quite  stiff,  with  bulk 
moduli  about  30%  smaller  than  diamond  while  the  shear 
modulus  is  about  60%  of  diamond.  This  is  somewhat  stiffer 
than  the  moduli  predicted  by  Aydin  et  al.  [30]  for  their  ReP4 
structure,  and  nearly  as  stiff  as  the  N2W  structures  predicted  by 
Wang  et  al.  [26] .  While  this  is  no  guarantee  that  the  structure  is 
in  fact  hard,  it  does  indicate  that  the  NbO  structure  of  tungsten 
nitride  is  a  candidate  for  a  superhard  material.  Unlike  the  N2W 
and  N4W  candidate  structures  mentioned  elsewhere,  every 


184110-7 


MEHL,  FINKENSTADT,  DANE,  HART,  AND  CURTAROLO 

TABLE  III.  Equilibrium  lattice  and  elastic  constants  of  tungsten 
nitride  in  the  NbO  structure  [space  group  Pm3m  #221,  Wyckoff 
positions  (3c)  and  (3d),  #53  in  the  Supplemental  Material  [48]]. 
These  were  computed  by  VASP  using  the  appropriate  PAW  potentials 
for  each  exchange-correlation  functional.  Elastic  constants  (in  GPa) 
were  computed  by  finite  strain  [80,81].  The  isotropic  shear  modulus 
G  is  the  average  of  the  Hashin-Shtrikman  bounds  for  a  cubic  system 
[88,89].  Starred  results  are  from  the  paper  of  Liu,  Zhou,  Gall,  and 
Khare  [67]. 
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functional  predicts  the  NbO  structure  to  be  the  ground- state 
structure  at  this  composition. 

In  addition  to  the  elastic  constants  shown  in  Table  III,  we 
further  checked  the  stability  of  the  NbO  structure  by  computing 
the  phonon  spectra,  as  shown  in  Fig.  5.  There  are  no  imaginary 
frequency  modes,  indicating  that  the  structure  is  stable. 

The  electronic  band  structure  and  density  of  states  of  the 
NbO  phase  are  shown  in  Figs.  6  and  7,  respectively.  As  in  all 
cases  in  this  paper,  the  high- symmetry  k  points  are  labeled 
following  as  in  Lax  [91]  and  by  AFLOW  [92].  The  structure  is 
a  metal,  with  most  of  the  electrons  at  the  Fermi  level  being 
provided  by  the  tungsten  d  states. 

The  change  in  enthalpy  AH  for  each  of  our  -phase 
structures  is  shown  in  Fig.  8.  This  is  a  subset  of  Fig.  2 
containing  only  /3 -phase  structures.  Each  point  is  color  coded 
to  show  the  concentration  of  vacancy  sites  on  the  underlying 
NaCl  lattice,  e.g.,  the  NbO  structure  (#53)  is  a  supercell  of 


r  x  r  mr  r 

FIG.  5.  (Color  online)  Phonon  frequencies  along  high- symmetry 
lines  for  the  NbO  structure  (#53)  of  tungsten  nitride,  found  with 
the  MedeA®  package  [79].  No  imaginary  phonon  frequencies  were 
found,  indicating  that  this  structure  is  stable  against  small  amplitude 
vibrations.  The  high- symmetry  points  and  lines  are  described  in  the 
Supplemental  Material  [48]. 
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FIG.  6.  (Color  online)  The  electronic  band  structure  of  tungsten 
nitride  in  the  NbO  structure  (#53)  at  the  PBE  equilibrium  ( a  = 
4.131  A).  The  high-symmetry  points  and  lines  are  described  in  the 
Supplemental  Material  [48]. 

the  NaCl  structure  with  25%  of  its  sites  vacant.  We  see  that 
the  lower-energy  structures,  particularly  in  the  region  of  50% 
tungsten,  are  dominated  by  unit  cells  with  25%  or  fewer 
vacancies.  This  suggests  that  the  experimentally  observed 
/LNW  structures  have  the  atoms  arranged  in  a  sodium  chloride 
structure  where  each  site  is  approximately  75%  occupied, 
the  exact  vacancy  fraction  for  each  site  depending  on  the 
composition.  This  picture  fits  with  the  experimental  lattice 
constant  [22]  for  the  the  /3  phase,  which  at  4.12-4.14  A  is 
similar  to  the  lattice  constants  we  obtain  for  the  NbO  structure 
(4.08  A  LDA,  4.13  A  PBE),  and  far  below  the  lattice  constant 
we  obtain  for  the  fully  occupied  NaCl  structure  (4.30  and 
4.37  A,  respectively). 

Contrary  to  published  results  [14-16,18,19,23],  we  see 
no  evidence  of  stable  /3 -phase  NW2  structures,  nor  of  any 
structures  where  all  of  the  tungsten  sites  are  occupied.  For  x 


e-eF  (eV) 

FIG.  7.  (Color  online)  The  electronic  density  of  states  and  angu¬ 
lar  momentum  decomposed  density  of  states  for  tungsten  nitride  in 
the  NbO  structure  (#53).  These  were  computed  via  the  tetrahedron 
method  in  VASP.  The  primary  contributions  are  from  the  N-p  and  W -d 
states,  with  the  W- d  states  dominating  the  conduction  band. 
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FIG.  8.  (Color  online)  A  section  of  Fig.  2  showing  only  structures  created  by  removing  atoms  from  supercells  of  the  NaCl  structure.  Yellow 
circles  denote  structures  where  25%  of  the  sites  on  the  NaCl  lattice  are  empty.  Downward-pointing  triangles  have  fewer  than  25%  of  the  sites 
vacant,  and  the  redder  triangles  have  fewer  vacancies,  culminating  in  the  NaCl  structure  (#67).  Upward-pointing  triangles  have  more  than  25% 
of  the  sites  vacant,  and  are  bluer  as  more  sites  become  vacant.  Known  experimental  structures  are  labeled  by  their  prototype  name. 


near  |,  the  closest  structure  -phase  structure  to  the  tie  line  is 
a  NgWi5  structure  (#88),  which  has  a  vacancy  on  the  tungsten 
site  as  well  as  multiple  vacancies  on  the  nitrogen  sites  of  the 
NaCl  supercell.  This  suggests  that  vacancies  on  both  nitrogen 
and  tungsten  sites  are  important  in  stabilizing  any  cubic  (/ 3 ) 
tungsten  nitride  structure. 

C.  8 -phase  structures 

While  the  (and  y)  phases  in  the  N-W  system  can  all  be 
described  by  various  patterns  of  both  nitrogen  and  tungsten 
vacancies  in  supercells  of  the  sodium  chloride  structure,  8 
phases  consist  of  alternating  triangular  sheets  of  nitrogen 
and  tungsten  atoms,  stacked  in  hexagonal  or  rhombohedral 
patterns,  as  in  tungsten  carbide  and  nickel  arsenide,  with 
experiments  [18,19,21]  suggesting  that  vacancies  appear  only 
on  the  tungsten  sublattices. 

Schonberg  [49]  suggested  that  the  parent  8  phase  is  the 
tungsten  carbide  structure  (#61  in  the  Supplemental  Material 
[48]),  with  vacancies  on  the  nitrogen  sites.  He  did  not  pursue 
this  idea,  as  at  the  time  the  determination  of  nitrogen  positions 
in  a  tungsten  nitride  compound  was  impossible. 

The  first  large-scale  investigation  of  8  phases  was  done  by 
Khitrova  and  Pinsker  [18].  They  found  six  5-phase  structures, 
ranging  in  composition  from  N2W  (5^,  structure  #25  in  the 
Supplemental  Material  [48])  to  NW2  (5^7,  #101).  While  these 
end  points  do  not  have  any  vacancies,  the  intermediate  struc¬ 
tures  (5^///,/y,  5^7)  have  vacancies  on  some  tungsten  planes. 

Recently,  WYL+  [21]  found  experimental  evidence  for  a  8- 
type  rhombohedral  structure  which  they  designated  ‘V-W2N3” 


(x  =  0.4).  This  structure  consists  of  close-packed  planes 
of  alternating  layers  of  tungsten  and  nitrogen,  with  the 
stacking  pattern  ABABCACABCBC.  Alternating  planes  of 
tungsten  have  an  occupation  factor  of  |,  i.e.,  two-thirds  of 
the  tungsten  sites  are  vacant,  giving  the  N3W2  stoichiometry. 
Interestingly,  WYL+  also  considered  a  structure  which  they 
called  “/1-W2N3”  (x  =  0.6,  #80);  this  structure  is  identical  to 
the  parent  (i.e.,  no  tungsten  vacancies)  8!H  structure  of  Khitrova 
and  Pinsker  (#80),  except  that  placement  of  the  tungsten  and 
nitrogen  atoms  is  reversed. 

As  we  saw  with  the  /3 -phase  structures,  atomistic  calcu¬ 
lations  cannot  readily  model  randomly  placed  vacancies.  We 
approximated  the  random  structures  by  creating  supercells  of 
the  parent  structures,  and  then  removing  various  patterns  of 
tungsten  atoms  to  match  the  required  stoichiometry.  We  also 
looked  at  structures  slightly  off  the  target  stoichiometry.  The 
Supplemental  Material  [48]  describes  the  low-lying  (within 
0. 2-0.3  meV /atom  of  the  convex  hull)  structures  we  found  by 
this  method.  It  is  necessarily  incomplete,  but  we  believe  we 
have  found  structures  extremely  close  to  the  actual  minimum- 
energy  structures. 

Figure  9  shows  some  of  the  low-energy  structures  in  this 
system,  as  a  function  of  tungsten  concentration  (along  the  v 
axis)  and  vacancy  concentration  (shading).  We  started  by  fully 
occupying  all  the  tungsten  sites  in  the  six  8  phases  described 
by  Khitrova  and  Pinsker  [18,19],  the  “r-W2N3”  phase  found 
by  WYL+  [21]  (the  circled  structures  in  the  diagram),  and 
the  WC  structure  (pentagons)  proposed  by  Schonberg  [49]. 
We  then  systematically  removed  atoms  from  the  appropriate 
tungsten  layers.  In  most  cases,  this  lowered  A H,  and  two 
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FIG.  9.  (Color  online)  Enthalpy  per  atom  for  various  crystal  structures  approximating  the  8  phases  of  Khitrova  and  Pinkser  [18,19]  and 
Wang,  Yu,  Lin  et  al.  (WYL+)  [21].  All  calculations  were  done  using  the  PBE  functional.  The  legend  describes  the  six  structures  found  by 
Khitrova  and  Pinsker,  as  well  as  the  ‘V-W2N3”  structure  of  WYL+.  We  also  include  the  tungsten  carbide  (WC)  structure  and  some  derivatives 
[25,49].  The  actual  structures  are  labeled  according  to  the  indexing  in  the  Supplemental  Material  [48].  The  circled  structures  are  the  parent 
structures,  which  have  fully  occupied  tungsten  sites.  For  the  other  structures,  the  tungsten  vacancy  concentration  is  indicated  by  the  point 
shading.  Red  (dark)  shaded  structures  have  no  or  few  defects,  while  orange-yellow  (light)  shading  indicates  a  larger  concentration  of  vacancies, 
up  to  50%.  Note  that  some  tungsten  planes  have  no  vacancies,  while  other  planes  can  have  as  many  as  |  of  the  sites  vacant.  See  Fig.  10  for  an 
example.  The  8!H  and  8 ^  structures  have  no  vacancy  sites,  therefore,  the  parent  is  the  only  structure  we  studied  and  its  special  status  is  shown 
by  having  no  coloring  to  the  data  point. 


of  our  derivative  ‘V-W2N3”  structures  (#35  and  #36)  have 
the  lowest  AH  of  all  structures  in  this  study,  independent  of 
choice  of  functional.  In  addition,  the  two  structures  (which  are 
numerically  degenerate  in  energy)  form  part  of  the  convex  hull 
for  the  N-W  system.  The  “/z-W2N3”  phase  mentioned  above 
is  approximately  0.35  eV/atom  above  the  ground-state  phase, 
so  it  is  at  best  a  metastable  state. 

Of  the  two  ground- state  structures,  #35  has  a  smaller  unit 
cell  than  #36,  so  we  will  concentrate  on  it.  Figure  10  shows  the 
atomic  positions  for  this  structure.  Its  elastic  constants,  listed  in 
Table  IV,  were  computed  using  the  finite  displacement  method 
contained  in  VASP,  and  were  used  to  calculate  the  Reuss  (lower) 
and  Voigt  (upper)  bounds  on  the  isotropic  poly  crystalline  bulk 
and  shear  moduli  [55] .  The  shear  modulus  is  less  than  200  GPa, 
so  we  cannot  regard  this  as  a  candidate  for  a  superhard 
structure. 

The  electronic  band  structure  is  shown  in  Fig.  11,  and  the 
electronic  density  of  states  in  Fig.  12.  There  is  a  small  density 
of  states  at  the  Fermi  level,  so  this  structure  is  a  semimetal. 
We  only  show  the  results  for  the  PBE  calculation,  but  this 
behavior  persists  across  all  three  choices  of  DFT.  Similar 
behavior  occurs  in  the  nearly  degenerate  structure  #36. 

Another  low-energy  8 -type  structure  is  the  hP  3  N2W 
structure  proposed  by  WLLX+  [26],  which  can  be  constructed 
by  doubling  the  tungsten  carbide  (#61)  unit  cell  in  the  z 


direction  and  removing  one  of  the  tungsten  atoms.  This 
structure  is  composed  of  N2  dimers  alternating  with  layers  of 
tungsten,  and  so  we  will  discuss  it  in  the  next  section,  although 
we  do  include  it  in  Fig.  9. 

In  addition  to  the  structures  described  above,  we  also  looked 
at  the  structure  labeled  “/j-W2N3”  by  WYL+  [21],  which  is 
same  as  the  8 !H  N2W3  structure  #80  with  N  and  W  occupations 
reversed,  forming  structure  #43.  The  aSm  structure  [93], 
structures  #25  and  #91,  consists  of  alternating  close-packed 
layers  (stacking  NNW  or  NWW),  and  so  can  also  be  considered 
a  8  phase.  Finally,  we  tried  modeling  5-like  structures  starting 
with  N-W  in  the  NaCl  structure,  considering  it  as  a  close- 
packed  system  with  stacking  ABCABCABC,  removing  two- 
thirds  of  the  atoms  in  alternating  tungsten  layers,  forming 
structure  #38.  While  these  structures  had  AH  <  0,  they  were 
not  near  the  convex  hull  using  any  functional. 

D.  N2  phases 

For  our  purposes,  an  N2  phase  is  one  in  which  we  can  readily 
identify  nitrogen  dimers  inside  a  tungsten  matrix.  Ay  din, 
Ciftci,  and  Tatar  (ACT)  [30]  studied  N4 W  in  the  ReP4  structure 
[63]  (Supplemental  Material  structure  #7)  They  equilibrated 
the  system  using  the  LDA  functional  and  found  a  negative 
value  for  AH,  and  also  computed  the  elastic  constants  and  an 
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FIG.  10.  (Color  online)  The  PBE  minimum-energy  structure 
#35,  space  group  Cm  —  C3X  (#8).  This  is  an  approximation  to  the 
“r-W2N3”  structure,  and  has  the  proper  number  of  tungsten  vacancies 
in  the  unit  cell  [21].  In  the  experimental  ‘V-W2N3”  structure,  the 
tungsten  sites  in  the  z  =  0  plane  are  randomly  occupied  and  do  not 
necessarily  have  the  pattern  shown  here. 

estimate  of  the  hardness,  stating  that  this  material  should  be 
superhard. 

In  the  form  proposed  by  ACT,  this  structure  is  of  type 
N2,  with  elongated  nitrogen  dimers  in  a  tungsten  matrix. 
It  is  extremely  difficult  to  equilibrate.  Starting  from  the 
structure  which  the  authors  graciously  provided  us,  we  found 
a  minimum-energy  structure  (#7)  close  to  theirs  as  shown  in 
Fig.  1.  Unfortunately,  we  found  that  AH  >  0,  suggesting  a 
difference  in  end-point  energies  in  Eq.  (1)  in  our  calculation 
compared  to  ACT.  In  addition,  when  we  used  the  ACT  structure 


FIG.  1 1 .  (Color  online)  The  electronic  band  structure  for  struc¬ 
ture  #35,  one  of  our  approximations  to  the  “r-W2N3”  structure  [21], 
using  the  PBE  functional. 
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TABLE  IV.  Density  functional  dependent  elastic  constants  of 
structure  #35,  space  group  Cm  (#8),  Pearson  symbol  mC 20,  shown 
in  Fig.  10.  This  is  an  approximate  version  of  the  ‘V-W2N3”  structure 
[21].  We  also  include  the  Reuss  (lower)  and  Voigt  (upper)  bounds  on 
the  isotropic  polycrystalline  bulk  and  shear  moduli  for  these  systems 
[55].  All  distances  are  in  A,  while  the  moduli  are  in  GPa.  Unlisted 
elastic  constants  are  zero  by  symmetry.  Note  that  the  Reuss  bulk 
modulus  Br  is  the  equilibrium  bulk  modulus  of  the  crystal  given  by 
K0  =  —VqP'(Vo),  where  Vo  is  the  equilibrium  volume  of  the  crystal. 
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FIG.  12.  (Color  online)  The  electronic  density  of  states  for  struc¬ 
ture  #35,  one  of  our  approximations  to  the  “r-W2N3”  structure  [21], 
using  the  PBE  functional  and  the  tetrahedron  method.  Note  that  the 
angular  momentum  decomposed  densities  of  states  are  summed  over 
all  atoms  of  the  same  type,  even  if  they  have  different  crystallographic 
positions. 
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with  the  PBE  and  vdW-DF2  functionals,  we  found  AH  to  be 
large  and  positive,  over  0.7  eV/atom  in  the  PBE  case. 

Furthermore,  with  some  manipulation  of  the  structures  we 
were  able  to  obtain  even  lower-energy  structures  (#6).  In  the 
LDA  case,  we  swapped  the  positions  of  the  tungsten  atoms 
with  the  second  set  of  nitrogen  atoms,  and  found  an  A  H  <0 
and  much  closer  to  the  convex  hull,  though  still  not  on  it.  This 
procedure  did  not  lead  to  a  low-energy  structure  for  the  PBE  or 
vdW-DF2  functions,  but  when  we  expanded  the  volume  of  the 
LDA  structure  we  found  a  double-well  structure  in  the  energy, 
where  the  higher  volume  had  a  lower  energy.  Backtracking 
from  this,  we  found  an  even  lower-energy  structure  at  a  volume 
comparable  to  the  original  LDA/ACT  structure.  When  we 
used  this  minimized  structure  with  the  PBE  and  vdW-DF2 
functionals,  we  obtained  the  structures  designated  (#6)  in  the 
Supplemental  Material  [48]. 

Given  these  results,  we  must  conclude  that  the  ACT 
structure  is  not  the  lowest-energy  structure  related  to  ReP4, 
that  is,  with  space  group  Pbca  and  all  atoms  on  (8c)  Wyckoff 
sites.  In  addition,  all  of  the  ReP4  structures  are  least  0.1  eV 
above  the  convex  hull  of  the  WN  system,  and  so  unlikely  to 
be  seen.  However,  we  are  not  certain  that  we  have  found  the 
minimum-energy  structure  within  these  constraints,  for  any 
choice  of  density  functional,  and  we  are  continuing  the  search. 

In  addition  to  the  ReP4  structure,  we  looked  other  candi¬ 
dates  for  the  N4W  ground  state.  The  only  low-energy  structure 
we  found  was  related  to  the  metastable  FeB4  structure  pro¬ 
posed  by  Van  der  Geest  and  Kolmogorov  [70].  We  originally 
found  the  enthalpy  of  this  structure  to  have  AH  ~  1  eV 
using  the  PBE  functional.  On  transferring  this  structure  from 
one  computer  to  another  we  inadvertently  transposed  some 
atomic  positions.  The  resulting  structure,  #8,  is  much  lower  in 
energy  than  the  original  converged  structure,  independent  of 
the  choice  of  functional,  but  it  is  still  far  from  the  convex  hull. 

Given  that  we  have  found  two  structures  with  relatively  low 
enthalpies,  we  cannot  eliminate  the  possibility  that  a  stable 
N4W  compound  exists,  but  at  the  moment  this  seems  unlikely. 

Compounds  with  stoichiometry  N2W  have  been  studied  by 
other  workers.  WLLX+  [26]  used  an  evolutionary  method  and 
found  two  nearly  degenerate  structures  which  we  designate 
as  #11  and  #12.  These  are  N2  dimers  stacked  in  a  hexagonal 
crystal  alternating  with  tungsten  atoms.  The  structures  differ  by 
the  stacking  patterns  of  the  dimers,  and  they  could  be  classified 
as  either  N2  phases  or  as  8  phases  with  missing  tungsten  planes. 

Table  V  gives  the  equilibrium  lattice  constants  and  elastic 
constants  of  the  two  structures  using  all  three  density  func¬ 
tionals.  We  compare  our  results  to  those  of  WLLX-I-,  who 
looked  at  the  systems  using  LDA  and  PBE.  In  all  cases,  both 
structures  have  relatively  large  bulk  and  shear  moduli,  and  so 
are  candidate  superhard  materials. 

The  structures  #11  and  #12  are  the  lowest-energy  N2W 
systems  within  our  LDA  and  PBE  calculations.  In  the  LDA 
they  form  part  of  the  convex  hull  for  the  tungsten  nitrogen 
system,  as  shown  in  Fig.  1. 

E.  Si02  phases 

One  of  our  approximate  /3 -phase  constructs,  which  we  will 
designate  cl 36  (#15  in  the  Supplemental  Material  [48]),  had 
a  surprisingly  low  enthalpy  when  we  calculated  its  energy 
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TABLE  V.  Equilibrium  structural  parameters  for  the  hexagonal 
N2W  structures  found  by  Wang  et  al.  (hereafter  WLLX+)  [26].  The 
low- symmetry  hP 3  structure  (#1 1  in  the  Supplemental  Material  [48]) 
is  in  space  group  P6m2,  with  nitrogen  atoms  at  the  (2 g)  Wyckoff 
positions  (00z)  and  the  tungsten  atom  on  the  (id)  site  (||  |).  In  the 
high-symmetry  hP 6  structure  (#12,  P63/mmc),  the  nitrogen  atoms 
are  on  ( 4e )  sites  (00z),  and  the  tungsten  atoms  are  on  (2d)  sites  (|  1 1). 
For  the  LDA  and  PBE  functionals,  we  show  our  results  (O),  along 
with  the  results  of  WLLX+.  Rn-n  and  RN.W  are  the  length  of  the 
respective  bonds.  The  elastic  constants  are  in  GPa.  The  isotropic  bulk 
(B)  and  shear  (G)  moduli  are  averages  of  the  Reuss  and  Voigt  bounds 
[55]. 


Functional 

LDA 

Low- symmetry 

PBE 

P6m2,  hP3 

vdW-DF2 

O 

O 

W 

O 

W 

a  (A) 

2.887 

2.887 

2.928 

2.933 

3.000 

c(A) 

3.877 

3.877 

3.916 

3.918 

3.974 

z 

0.1804 

0.1804 

0.1813 

0.1814 

0.1824 

Rn-n  (A) 

1.399 

1.399 

1.421 

1.421 

1.450 

Rn-w  (A) 

2.077 

2.077 

2.104 

2.104 

2.143 

G11 

638 

654 

573 

588 

493 

C33 

1056 

1082 

954 

973 

827 

C44 

241 

260 

222 

232 

191 

C\2 

207 

213 

183 

191 

159 

C13 

230 

248 

193 

206 

157 

B 

396 

412 

351 

255 

297 

G 

246 

255 

226 

231 

195 

High- symmetry  P  63 /mmc,  hP6 

O 

W 

O 

W 

O 

a  (A) 

2.893 

2.893 

2.939 

2.939 

3.007 

c(A) 

7.714 

7.714 

7.796 

7.796 

7.891 

z 

0.0898 

0.0898 

0.0902 

0.0902 

0.0907 

RN-n  (A) 

1.392 

1.392 

1.406 

1.406 

1.431 

RN-w  (A) 

2.078 

2.078 

2.105 

2.105 

2.144 

C11 

633 

642 

568 

579 

486 

C33 

1051 

1078 

952 

973 

827 

C44 

249 

262 

228 

233 

197 

c  12 

213 

217 

187 

195 

161 

C13 

237 

252 

199 

211 

197 

B 

399 

411 

353 

364 

298 

G 

245 

252 

225 

228 

195 

using  the  van  der  Waals  vdW-DF2  functional.  The  cl 36 
structure  starts  from  a  32-atom  body-centered-cubic  supercell 
of  the  NaCl  structure,  with  4  nitrogen  atoms  and  10  tungsten 
atoms  removed.  This  leaves  only  the  (24 h)  and  (12 d)  Wyckoff 
positions  occupied,  and  so  we  originally  classified  it  as  a 
ft -phase  structure.  However,  as  seen  in  Fig.  13  and  described  in 
Table  VI,  each  tungsten  atom  is  bonded  to  four  nitrogen  atoms 
arranged  tetrahedrally,  and  each  nitrogen  atom  is  bound  to  two 
tungsten  atoms.  This  arrangement  is  reminiscent  of  many  SiC>2 
structures  [69],  although  this  particular  structure  has  not  been 
seen  experimentally. 

We  also  looked  some  of  the  simpler  Si02  structures, 
including  the  high-temperature  phase  of  quartz  [94]  (#16  in 
the  Supplemental  Material  [48]),  cristobalite  [95]  (#14),  and 
tridymite  [96]  (#13).  The  low-temperature  forms  of  these 
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TABLE  VI.  The  cl 36  structure  of  N2W,  constructed  from  a  32- 
atom  body-centered-cubic  supercell  of  the  NaCl  structure.  The  space 
group  is  Im3m-Ol  (#229),  the  nitrogen  atoms  occupy  the  (24 h) 
Wyckoff  positions  (Oyy),  and  the  tungsten  atoms  occupy  the  (12 d) 
sites  (0").  The  first  N-W-N  angle  is  for  atoms  lying  on  the  surface  of 
the  cube  cell.  The  second  case  has  one  of  the  N  atoms  in  the  interior 
of  the  cube. 


Functional 

LDA 

PBE 

vdW 

a  (A) 

10.273 

10.383 

10.490 

y 

0.1454 

0.1455 

0.1455 

N-W  bond  (A) 

1.840 

1.860 

1.879 

Angles 

W-N-W 

108.7° 

108.6° 

108.5° 

N-W-N  (face) 

108.7° 

108.6° 

108.5° 

N-W-N  (interior) 

109.9° 

109.9° 

110.0° 

structures  are  degenerate  with  these,  so  we  do  not  report  them 
here.  At  x  =  |  the  high- temperature  quartz  structure  had  the 
lowest  enthalpy  when  using  the  van  der  Waals-DF2  functional, 
and  it  forms  part  of  the  ground-state  hull  in  Fig.  3. 

We  also  looked  at  some  distorted  rutilelike  structures,  which 
have  bonding  similar  to  SiC>2.  These  include  bixbyite  [97,98], 
structures  #41  and  #42;  /INbC^  [99],  structure  #17;  and  PbC>2 
[100]  structure  #21,  which  has  the  same  valence  electron  count 
at  NbC>2.  These  structures  had  energies  comparable  to  the  SiC>2 
structures  discussed  above.  In  particular,  the  relaxed  /3NbC>2 
structure  is  quite  close  to  the  convex  hull  when  we  used  the 
vdW-DF2  functional. 

All  of  the  these  structures  are  insulating.  Since  the  cubic 
cl 36  structure  was  the  first  of  these  we  found,  and  has 
the  highest- symmetry  primitive  cell,  we  did  band- structure 
calculations  for  it,  as  shown  in  Fig.  14.  There  is  a  direct 


FIG.  13.  (Color  online)  The  Si02-like  cl 36  structure  (#15  in  the 
Supplemental  Material  [48])  predicted  as  the  ground  state  of  N2W  by 
the  vdW-DF2  functional.  The  basic  structure  looks  the  same  using 
LDA,  PBE,  and  vdW-DF2  functionals.  The  exact  dimensions  are 
described  in  Table  VI.  Note  that  this  is  a  body-centered-cubic  lattice, 
so  the  center  of  the  cube  and  the  cube  corners  are  equivalent  sites. 
The  inversion  sites  are  at  the  centers  of  the  eight-atom  rings,  not  on 
the  tungsten  atoms. 
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FIG.  14.  (Color  online)  The  electronic  band  structure  for  the 
Si02-like  body-centered-cubic  cl 36  structure  (#15)  of  N2W  de¬ 
scribed  in  the  text.  Calculations  were  done  using  the  vdW-DF2 
functional,  including  van  der  Waals  forces,  within  VASP.  The  top 
valence  band  has  a  very  small  dispersion,  dropping  to  —0.05  eV 
below  the  top  of  the  valence  band  at  P.  The  gap  is  0.75  eV,  and  is 
direct  at  T.  The  LDA  and  PBE  calculations  produce  similar  band 
structures. 

band  gap  of  0.75  eV  using  the  vdW-DF2  functional.  The  band 
structure  is  unusual  in  that  there  is  a  very  flat  band  (0.05  eV 
maximum  dispersion)  at  the  top  of  the  valence  band.  As  shown 
in  the  electronic  density  of  states  plot,  Fig.  15,  this  narrow  band 
has  strong  nitrogen  2 p  character. 

We  also  computed  the  elastic  constants  for  the  cl  36 
structure.  Not  surprisingly,  this  very  low-density  system  has 
small  elastic  constants,  as  shown  in  Table  VII.  We  expect  the 
other  SiC>2  structures  to  have  similar  shear  moduli. 

Given  the  large  density  of  states  near  the  top  of  the 
occupied  bands  in  this  system,  we  looked  for  structural 
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FIG.  15.  (Color  online)  The  electronic  total,  and  angular  momen¬ 
tum  decomposed,  density  of  states  for  N2W  in  the  Si02-like  cl 36 
structure  (#15),  using  the  vdW-DF2  functional  and  the  tetrahedron 
method  within  VASP.  We  truncated  the  plot  at  75  states/eV/cell 
in  order  to  show  the  behavior  of  the  system  near  the  top  of  the 
valence  band.  The  density  of  states  at  the  bottom  of  the  valence  band 
approaches  200  states/e V/cell. 
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TABLE  VII.  Elastic  constants  for  the  cl 36  structure  of  N2W, 
computed  using  the  LDA,  PBE,  and  vdW-DF2  functionals.  Elastic 
constants  were  computed  by  the  finite  strain  method  [80],  allowing 
the  atoms  to  relax  at  each  strain  while  fixing  the  unit  cell.  Equilibrium 
structural  parameters  are  taken  from  Table  VI.  B  is  the  equilibrium 
bulk  modulus,  and  all  elastic  constants  are  in  GPa.  The  shear  modulus 
G  is  an  average  of  the  Hashin-Shtrikman  bounds  [88,89]. 


Functional 

Cn 

C\2 

C44 

B 

Cn  —  C\2 

G 

LDA 

119 

77 

13 

91 

42 

15 

PBE 

112 

72 

13 

85 

41 

16 

vdW-DF2 

104 

65 

12 

78 

39 

14 

instabilities,  including  imaginary  phonon  frequencies  and  a 
possible  magnetic  state.  While  we  have  not  done  a  complete 
examination  of  the  phonon  spectrum,  those  k  points  we 
sampled  all  have  positive  phonon  frequencies,  and  imposing  a 
spin  polarization  on  the  crystal  raises  its  energy.  We  conclude 
that  the  structure  is  at  least  metastable. 

While  all  of  SiC>2  structures  we  examined  always  have 
AH  <  0,  none  of  them  form  part  of  the  convex  hull  when 
we  use  the  PBE  or  LDA  functionals.  The  appearance  of  these 
structures  on  and  near  the  vdW-DF2  convex  hull  (Fig.  3)  is 
most  likely  an  artifact  of  the  van  der  Waals  density  functional. 

F.  Miscellaneous  structures 

In  this  section,  we  describe  several  low-energy  structures 
that  appear  in  Figs.  1-3.  Many  of  these  structures  have  been 
studied  by  others,  as  outlined  in  Sec.  II. 

While  we  did  calculations  for  a  number  of  structures 
with  composition  NW2,  including  and  8  phases,  the  most 
interesting  structures  we  found  were  NW2  in  the  hexagonal 
M0S2  structure  [101]  (#89)  and  the  rhombohedral  q'MoS2 
structure  [102]  (#90).  In  all  cases,  the  hexagonal  M0S2 
structure  had  the  lowest  enthalpy  of  any  NW2  structure  we 
studied,  and  within  the  LDA  it  was  nearly  touching  the  tie 
line.  This  structure  has  not  been  seen  nor  previously  predicted 
in  the  N-W  system.  We  show  the  LDA  equilibrium  structure 
in  Fig.  16.  While  exhibiting  the  layered  character  of  the  M0S2 
structures,  it  is  not  two  dimensional.  In  fact,  the  distance 
between  atoms  in  adjacent  tungsten  layers  (2.81  A)  is  smaller 
than  the  distance  between  tungsten  atoms  in  the  plane  (2.84  A). 

This  structure  was  only  run  because  it  was  available  in 
the  AFLOW  database.  If  it  had  not  been  present  there,  it 
would  not  have  been  tested  since  we  had  no  indication  either 
experimentally  or  theoretically  that  the  N-W  system  would 
support  close-packed  tungsten  planes.  This  demonstrates  the 
power  of  the  structural  database  in  investigating  new  systems. 

Since  this  structure  has  a  similar  form  to  the  8  structures,  we 
tried  to  lower  its  energy  by  placing  vacancies  on  the  tungsten 
sites  (#85  in  the  Supplemental  Material  [48]  tables  and  in 
Figs.  1-3)  and  by  adding  interplaner  nitrogen  or  tungsten  atoms 
(#75,  #89).  In  no  case  were  we  able  to  get  an  enthalpy  below 
the  parent  M0S2  structure. 

This  structure  is  close  to  the  convex  hull  in  the  LDA,  and 
in  the  hopes  of  finding  an  instability  leading  to  a  lower-energy 
structure,  we  did  a  preliminary  test  for  hardness  by  computing 
the  elastic  constants  of  the  hexagonal  M0S2  phase,  shown  in 


FIG.  16.  (Color  online)  NW2  in  the  hexagonal  MoS2  structure 
[101]  #89  in  the  Supplemental  Material  [48].  Unlike  MoS2,  the 
interplaner  W-W  distance  is  smaller  than  the  in-plane  W-W  distance. 

Table  VIII.  It  is  not  particularly  hard,  having  a  shear  modulus 
on  the  order  of  100  GPa.  It  is,  however,  elastically  stable,  as  the 
elastic  constants  satisfy  the  Born  criterion  [39].  The  structure 
is  metallic,  as  can  be  seen  from  the  band  structure  (Fig.  17) 
and  electronic  density  of  states  (Fig.  18) . 

We  looked  at  various  structures  with  composition  NW, 
including  nickel  arsenide  (#59),  studied  by  Kroll,  Schroter, 

TABLE  VIII.  Density  functional  dependent  elastic  constants  of 
structure  of  NW2  in  the  hexagonal  MoS2  structure  [101]  #89  in  the 
Supplemental  Material  [48].  We  also  include  the  Reuss  (lower)  and 
Voigt  (upper)  bounds  on  the  isotropic  poly  crystalline  bulk  and  shear 
moduli  for  these  systems  [55] .  All  distances  are  in  A,  while  the  moduli 
are  in  GPa.  z  is  the  free  coordinate  of  the  tungsten  atom  at  the  (4/) 
Wyckoff  site.  Since  this  is  a  hexagonal  crystal,  C66  =  (C n  —  C i2)/2. 
Note  that  the  Reuss  bulk  modulus  BR  is  the  equilibrium  bulk  modulus 
of  the  crystal  given  by  K0  =  —VqP'(Vq). 


Functional 

Lattice  parameters 

LDA 

PBE 

vdW-DF2 

a  (A) 

2.837 

2.875 

2.933 

c  (A) 

10.212 

10.369 

10.584 

z[W  (4/)] 

0.612 

0.612 

0.611 

Elastic  constants  (GPa) 
Cn 

568 

512 

396 

C33 

730 

668 

581 

C44 

158 

142 

113 

C\2 

349 

312 

307 

C\3 

238 

207 

178 

Bulk  modulus  (GPa) 

Reuss  ( Br ) 

390 

349 

299 

Voigt  {By) 

391 

349 

300 

Shear  modulus  (GPa) 
Reuss  ( Gr ) 

142 

130 

74 

Voigt  (Gy) 

154 

141 

102 
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FIG.  17.  (Color  online)  The  electronic  band  structure  for  NW2 
in  the  hexagonal  MoS2  structure  [101,102]  (#89  in  the  Supplemental 
Material  [48]),  using  the  PBE  functional. 

and  Peters  [24],  and  the  zinc-blende  structure  (#65),  studied 
by  WYL+  [21].  The  nickel  arsenide  structure,  which  could  be 
considered  a  8  phase,  is  a  few  tenths  of  an  eV  above  the  NbO 
structure.  The  zinc-blende  structure,  which  WYL+  find  to  be 
very  stiff,  is  over  0.5  eV  above  NbO  and  so  is  unlikely  to  form. 

No  N3W  structures  were  found  near  the  convex  hull.  The 
closest  two  were  structures  molybdite  (MO3  [103],  #9)  and 
P3TC  (#10)  [56]  Song  and  Wang  [29]  studied  P3TC  within  the 
LDA,  and  found  to  have  reasonable  large  bulk  (303  GPa)  and 
shear  (156  GPa)  moduli.  They  also  found  that  it  had  negative 
enthalpy  (1),  which  we  confirm.  However,  as  shown  in  Fig.  1, 
this  structure  is  well  above  the  convex  hull  within  the  LDA. 
The  PBE  and  vdW-DF2  results  show  that  the  AH  for  P3TC  is 
positive.  It  is  unlikely,  therefore,  that  N3W  exists  in  the  P3TC 
structure. 

One-half  of  the  nitrogen  atoms  in  the  relaxed  molybdite 
structure  (which  looks  nothing  like  molybdite)  form  N2  dimers 
with  spacing  1.12-1.13  A.  The  relaxed  P3TC  structure  is 


e-eF  (eV) 


FIG.  18.  (Color  online)  The  electronic  total,  and  angular  momen¬ 
tum  decomposed,  density  of  states  for  NW2  in  the  hexagonal  MoS2 
structure  [101,102]  (#89  in  the  Supplemental  Material  [48]),  using 
the  PBE  functional. 


similar,  with  N2  dimers  spaced  1.2  A  apart.  We  can  therefore 
classify  both  structures  as  partial  N2  phases. 

Both  structures  have  the  same  Pnma  space  group,  with 
all  atoms  occupying  (8c)  Wyckoff  positions,  so  it  may  be 
possible  to  transform  from  P3TC  to  M0O3  by  simply  varying 
the  c/a  and  b /a  ratios.  Such  a  transformation  would  reveal  the 
structure  of  the  energy  barrier  between  the  two  states.  Since 
these  structures  are  well  above  the  minimum-energy  structure, 
we  have  not  investigated  this  further. 

V.  DISCUSSION 

To  truly  predict  the  structure  of  a  compound  AMB^v,  we 
must  find  the  lowest-energy  structure  at  that  composition, 
but  we  must  also  determine  if  that  structure  is  unstable  with 
respect  to  decomposition  into  either  into  its  end  points  or 
to  some  intermediate  structures.  When  doing  approximate 
density  functional  theory  calculations,  we  should  also  be 
sure  that  the  predicted  structure  is  stable  or  metastable  with 
respect  to  different  functionals  since  current  functionals  are 
only  approximations  to  the  real  one. 

Using  the  tungsten  nitride  system  as  an  example,  the  AFLOW 
high-throughput  method  was  used  to  calculate  the  relative 
enthalpy  (1)  of  over  500  structures  at  variety  of  compositions 
Ni_xWx.  We  then  studied  over  100  of  these  structures  in 
detail,  using  the  generalized-gradient  PBE  functional,  the 
local  density  approximation,  and  the  van  der  Waals  vdW-DF2 
functional  to  produce  a  zero-temperature  phase  diagram  over 
a  wide  variety  of  compositions. 

For  x^0.4,  all  three  functionals  give  the  same  predicted 
ground  states:  at  x  =  0.4,  a  hexagonal  8- NW  structure 
with  random  vacancies  on  the  tungsten  sites,  which  we 
approximated  using  various  supercells  of  the  parent  structure 
with  different  patterns  of  tungsten  vacancies;  and  the  /3-NW 
structure,  with  the  lowest  state  the  ordered  NbO  structure,  the 
sodium  chloride  structure  with  one-fourth  of  the  sites  vacant. 

All  three  functionals  show  the  importance  of  vacancies 
in  the  tungsten  nitride  system.  If  we  consider  the  sodium 
chloride  and  tungsten  carbide  structures  as  the  parents  of  the 
NW  /3  and  8  phases,  we  can  nearly  always  lower  the  energy 
by  constructing  a  supercell  and  removing  selected  atoms,  as 
shown  in  Figs.  8  and  9.  All  of  the  0  <  x  <  1  structures  on  the 
convex  hull  in  Figs.  1-3  can  be  constructed  by  removing  atoms 
from  supercells  of  the  parent  structures.  This  is  in  agreement 
with  experiment,  which  finds  random  vacancies  in  almost 
every  structure. 

What  is  not  in  agreement  with  experiment  is  the  actual 
composition  of  the  phase.  We  find  the  lowest-energy 
structures  near  x  =  ^ ,  while  experiment  tends  to  favor  cubic 
structures  with  x  ~  The  likeliest  explanation  is  in  the 
kinetics  of  the  tungsten-nitrogen  reaction.  It  is  very  difficult 
to  dissolve  N2  molecules  into  tungsten  [23],  and  so  it  is 
plausible  that  the  energy  barriers  in  the  problem  will  favor 
tungsten-rich  compounds.  A  thorough  investigation  of  the 
kinetics  of  tungsten  nitride  formation  is  necessary  for  further 
understanding. 

As  noted,  our  original  motivation  for  looking  at  the  tungsten 
nitride  system  was  the  possibility  of  forming  hard  carbon-free 
materials.  Previous  calculations  suggested  that  specific  N2W 
and  N4W  structures  were  stable,  had  large  bulk  and  shear 
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moduli  (over  200  GPa),  and  were  candidate  hard  structures. 
We  found  the  N4W  compound  [30]  to  be  well  above  the  convex 
hull,  so  that  it  would  tend  to  phase  separate  into  N2  and  either 
N2W  or  N3W2,  depending  on  the  choice  of  density  functional. 
The  candidate  N2W  structures  [26]  are  on  the  LDA  convex  hull, 
but  not  the  PBE  or  vdW-DF2  hulls,  and  so  might  be  found  at 
some  future  date.  The  experimentally  observable  ‘V-W2N3” 
structure  has  the  lowest  value  of  AH  in  the  entire  system. 
Unfortunately,  it  is  not  particularly  hard,  as  our  estimated  shear 
modulus  is  under  200  GPa.  However,  the  /3 -phase  ground- 
state  NbO  structure  does  have  shear  moduli  over  200  GPa, 
and  so  might  be  “hard.”  This  will  require  more  investigation, 
similar  to  what  Aydin,  Ciftci,  and  Tatar  [30]  did  for  their  N4W 
candidate. 

Another  surprise  coming  out  of  this  study  is  the  change 
in  the  shape,  depth,  and  composition  of  the  convex  hull 
depending  on  the  choice  of  density  functional.  In  the  absence 
of  spin-polarization  effects,  the  DFT  community  usually  thinks 
of  the  generalized-gradient  approximation,  and  especially 
the  PBE  and  related  functionals,  as  improving  on  the  EDA 
predictions  of  lattice  constants,  but  not  changing  the  ordering 
of  structures.  For  v  ^  0.4,  that  is  more  or  less  true  here. 
However,  for  nitrogen-rich  systems  there  is  considerable 
disagreement  about  the  ordering  of  structures.  This  is  most 
likely  caused  by  the  neglect  of  van  der  Waals  forces  in  the 
nitrogen-rich  systems.  To  test  this,  we  used  a  common  van 
der  Waals  functional,  vdW-DF2.  It  indeed  provides  a  better 
lattice  constant  for  the  solid  N2,  but  it  also  leads  to  what  seems 
to  be  a  spurious  prediction  of  the  ground- state  structure  at 
N2W.  It  also  consistently  overestimates  the  lattice  constants 
for  systems  with  jc  >  | ,  with  equilibrium  volumes  even  larger 
than  the  PBE  predictions. 

In  summary,  we  have  used  high-throughput  density  func¬ 
tional  calculations  to  map  out  the  low-energy  structures  in  the 
tungsten  nitride  system.  In  general,  we  find  that  previously 
proposed  structures,  both  theoretical  and  experimental,  are 
not  the  ground- state  structures  in  this  system,  with  the 
exception  of  the  “r-W2N3”  structure.  We  do  find  that  /Lphase 
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structures  can  be  described  as  various  vacancy  patterns  on 
the  sodium  chloride  structure,  but  we  cannot  explain  the 
experimental  preference  for  tungsten-rich  -phase  structures. 
We  also  find  that  none  of  the  density  functionals  we  used 
correctly  describe  the  system  over  all  compositions,  as  the 
LDA  and  PBE  neglect  van  der  Waals  forces  and  the  vdW-DF2 
functional  produces  what  appear  to  be  spurious  predictions  for 
equilibrium  structures. 

We  have  also  shown  that  predictions  of  structural  stability 
must  include  knowledge  of  structures  over  the  entire  range 
of  compositions.  Many  of  the  studies  described  in  Sec.  II 
found  structures  with  negative  formation  energy  AH  [Eq. 
(1)],  but  most  of  these  are  not  on  the  convex  hull  of  the 
tungsten  nitride  system.  One  of  the  structures  that  is  close 
to  the  hull,  M0S2,  would  never  have  been  studied  if  it  were 
not  in  the  database,  as  there  was  no  indication  that  it  was 
likely  to  occur.  The  high-throughput  mechanism  provided  by 
AFLOW  is  therefore  necessary  to  understand  any  compound 
system. 
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